##Cook, et al `Two Wrongs Make a Right'######
##Replication code for Illustration##########
##Last edit: 11/1/2016#######################

setwd("C:/Users/sjcook/Downloads/dataverse_files") #Specify working directory 
rm(list = ls())

install.packages("sp")
install.packages("maptools")
install.packages("dplyr")
install.packages("countrycode")
install.packages("cshapes")
install.packages("zoo")
install.packages("pracma")
install.packages("geosphere")
install.packages("matrixStats")
install.packages("dummies")
install.packages("reshape")
install.packages("data.table")
install.packages("lubridate")

library(sp)
library(maptools)
library(dplyr)
library(countrycode)
library(cshapes)
library(zoo)
library(pracma)
library(geosphere)
library(matrixStats)
library(dummies)
library(reshape)
library(data.table)
library(lubridate)

##Need old version of McSpatial##
require(devtools)
install_url("http://CRAN.R-project.org/src/contrib/Archive/McSpatial/McSpatial_1.1.1.tar.gz")
library(McSpatial)

##Construct Sample Dimensions##

source("state_panel.R")

state.panel <- state_panel("2012-01", "2013-12", by="month")

state.panel$parse <- parse_date_time(state.panel$date,"ymd") 

state.panel$styr <- year(state.panel$parse)
state.panel$stmo <- month(state.panel$parse)


##Read in/Arrange SCAD event data##
scad <-read.csv("SCAD.csv")
scad = filter(scad, etype == 7) #repression

scad.panel <- merge(state.panel, scad, by = c("ccode","styr","stmo"),  all.x = T)

scad.panel$cont<-countrycode(scad.panel$ccode,"cown","continent")
afr.scad.panel= filter(scad.panel, cont=="Africa")

afr.scad.panel$AP = ifelse(afr.scad.panel$nsource %in% c('AP','Both'),1,0)
afr.scad.panel$AFP = ifelse(afr.scad.panel$nsource %in% c('AFP','Both'),1,0)

afr.scad.panel.agg <- aggregate(x = cbind(afr.scad.panel$AP,afr.scad.panel$AFP), by = list(afr.scad.panel$ccode,afr.scad.panel$styr,afr.scad.panel$stmo), FUN = "sum")

names(afr.scad.panel.agg)[names(afr.scad.panel.agg) == 'Group.1'] <- 'ccode'
names(afr.scad.panel.agg)[names(afr.scad.panel.agg) == 'Group.2'] <- 'styr'
names(afr.scad.panel.agg)[names(afr.scad.panel.agg) == 'Group.3'] <- 'stmo'
names(afr.scad.panel.agg)[names(afr.scad.panel.agg) == 'V1'] <- 'AP'
names(afr.scad.panel.agg)[names(afr.scad.panel.agg) == 'V2'] <- 'AFP'

afr.scad.panel.agg$AP.dum <- ifelse(afr.scad.panel.agg$AP>0,1,0)
afr.scad.panel.agg$AFP.dum <- ifelse(afr.scad.panel.agg$AFP>0,1,0)

xtabs(~AP.dum+AFP.dum, data=afr.scad.panel.agg)

reports<-read.table("Africa.AP.AF.Reports.2012.2013.tab", header=T)
names(reports)[names(reports) == 'AP'] <- 'ap.reports'
names(reports)[names(reports) == 'AF'] <- 'afp.reports'
names(reports)[names(reports) == 'year'] <- 'styr'
names(reports)[names(reports) == 'month'] <- 'stmo'
names(reports)[names(reports) == 'CCode'] <- 'ccode'

rep.data <- merge(afr.scad.panel.agg, reports, by = c("ccode","styr","stmo"),  all.y = T)

rep.data$both.dum = as.numeric(rep.data$AP.dum+rep.data$AFP.dum>=1)

##Read in Predictor Data##
predictors <- read.table("predictors.tab", header=T, fill=T)

names(predictors)[names(predictors) == 'year'] <- 'styr'

est.data <- left_join(rep.data,predictors,by=c("ccode","styr"))

est.data <- merge(rep.data,predictors,by=c("ccode","styr"),  all.x = T)

est.data <- na.omit(est.data)

##Estimate models that are presented in Tables 3&4 (below)##

m1 <- summary(glm(est.data$both.dum ~ est.data$lag_gdppc_ln + est.data$lag_pop_ln + est.data$lag_demo  , family=binomial(link="probit")))

m2 <- misclass(est.data$both.dum ~ est.data$lag_gdppc_ln + est.data$lag_pop_ln + est.data$lag_demo , a0=FALSE) 

source("misclass_phi1_kpred.R")

m3 <- misclass_phi1_kpred(est.data$both.dum ~ est.data$lag_gdppc_ln  + est.data$lag_pop_ln + est.data$lag_demo  + est.data$ap.reports + est.data$afp.reports) 

source("misclass_a1_a2.R")

m4 <- misclass_a1_a2(est.data$both.dum ~ est.data$lag_gdppc_ln + est.data$lag_pop_ln + est.data$lag_demo , est.data$AP.dum, est.data$AFP.dum)

source("misclass_a1xz_a2xz.R")

m5 <-misclass_a1xz_a2xz(est.data$both.dum ~ est.data$lag_gdppc_ln + est.data$lag_pop_ln + est.data$lag_demo 
                   + est.data$ap.reports + est.data$afp.reports, est.data$AP.dum, est.data$AFP.dum)



##########Table 3#############

beta.gdp <- c(m1$coef[2],m2$estimate[2],m3$estimate[2], m4$estimate[2], m5[1,2])
se.gdp <- c(m1$coef[2,2],m2$stderr[2],m3$stderr[2], m4$stderr[2], m5[2,2] )

beta.pop <- c(m1$coef[3],m2$estimate[3],m3$estimate[3], m4$estimate[3], m5[1,3])
se.pop <- c(m1$coef[3,2],m2$stderr[3],m3$stderr[3], m4$stderr[3], m5[2,3] )

beta.demo <- c(m1$coef[4],m2$estimate[4],m3$estimate[4], m4$estimate[4], m5[1,4])
se.demo <- c(m1$coef[4,2],m2$stderr[4],m3$stderr[4], m4$stderr[4], m5[2,4] )

beta.cons <- c(m1$coef[1],m2$estimate[1],m3$estimate[1], m4$estimate[1], m5[1,1])
se.cons <- c(m1$coef[1,2],m2$stderr[1],m3$stderr[1], m4$stderr[1], m5[2,1] )

t3 <-matrix(c(beta.gdp,se.gdp,beta.pop,se.pop,beta.demo,se.demo,beta.cons,se.cons),ncol=5,byrow=T)
colnames(t3) <- c("M1","M2","M3","M4","M5")
rownames(t3) <- c("beta.gdp","se.gdp","beta.pop","se.pop","beta.demo","se.demo","beta.cons","se.cons")
t3 <- round(t3,3)
t3

##########Table 4#############

###Top Section###

beta.ap.gdp <- c(0,0,m3$estimate[6],0, m5[1,6])
se.ap.gdp <-  c(0,0,m3$stderr[6],0, m5[2,6])

beta.ap.pop <- c(0,0,m3$estimate[7],0, m5[1,7])
se.ap.pop <- c(0,0,m3$stderr[7],0, m5[2,7])

beta.ap.demo <- c(0,0,m3$estimate[8],0, m5[1,8])
se.ap.demo <- c(0,0,m3$stderr[8],0, m5[2,8])

beta.ap.reports <- c(0,0,m3$estimate[9],0, m5[1,9])
se.ap.reports <- c(0,0,m3$stderr[9],0, m5[2,9])

beta.ap.cons <-  c(0,m2[2],m3$estimate[5],m4$estimate[5], m5[1,5])
se.ap.cons <- c(0,0,m3$stderr[5],m4$stderr[5], m5[2,5])

t4.ap <-matrix(round(as.numeric(c(beta.ap.gdp,se.ap.gdp,beta.ap.pop,se.ap.pop,beta.ap.demo,se.ap.demo,beta.ap.reports,se.ap.reports,beta.ap.cons,se.ap.cons)),3),ncol=5,byrow=T)
colnames(t4.ap) <- c("M1","M2","M3","M4","M5")
rownames(t4.ap) <- c("beta.gdp","se.gdp","beta.pop","se.pop","beta.demo","se.demo","beta.reports","se.reports","beta.cons","se.cons")
t4.ap[t4.ap == 0] <- NA
t4.ap ##Note that the S.E. for the constant in Model 3 is 0.003 not NA; there was no way to efficently extract this information)

###Bottom Section###

beta.afp.gdp <- c(0,0,0,0, m5[1,11])
se.afp.gdp <-  c(0,0,0,0, m5[2,11])

beta.afp.pop <- c(0,0,0,0, m5[1,12])
se.afp.pop <- c(0,0,0,0, m5[2,12])

beta.afp.demo <- c(0,0,0,0, m5[1,13])
se.afp.demo <- c(0,0,0,0, m5[2,13])

beta.afp.reports <- c(0,0,m3$estimate[10],0, m5[1,14])
se.afp.reports <- c(0,0,m3$stderr[10],0, m5[2,14])

beta.afp.cons <-  c(0,0,0,m4$estimate[6], m5[1,10])
se.afp.cons <- c(0,0,0,m4$stderr[6], m5[2,10])

t4.afp <-matrix(round(as.numeric(c(beta.afp.gdp,se.afp.gdp,beta.afp.pop,se.afp.pop,beta.afp.demo,se.afp.demo,beta.afp.reports,se.afp.reports,beta.afp.cons,se.afp.cons)),3),ncol=5,byrow=T)
colnames(t4.afp) <- c("M1","M2","M3","M4","M5")
rownames(t4.afp) <- c("beta.gdp","se.gdp","beta.pop","se.pop","beta.demo","se.demo","beta.reports","se.reports","beta.cons","se.cons")
t4.afp[t4.afp == 0] <- NA
t4.afp



